{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Upper Limits with Python\n",
    "\n",
    "This sample analysis shows a way to determine an upper limit on the GeV emission from Swift J164449.3+573451 similar to what was done in [Burrows et al. (Nature 476, page 421)](http://www.nature.com/nature/journal/v476/n7361/full/nature10374.html).\n",
    "\n",
    "To compute the upper limit, we use the profile likelihood method. This entails scanning in values of the normalization parameter, fitting with respect to the other remaining free parameters, and plotting the change in log-likelihood as a function of flux.\n",
    "\n",
    "Assuming 2$\\cdot$log-likelihood behaves asymptotically as chi-square, a 90% confidence region will correspond to a change in log-likelihood of 2.71/2.\n",
    "\n",
    "Note that this change in log-likelihood corresponds to a two-sided confidence interval. Since we are interested in an upper-limit, this change in log-likelihood actually corresponds to a 95% CL upper-limit. See [Rolke et al. (2005)](https://arxiv.org/abs/physics/0403059) for more details.\n",
    "\n",
    "We will first cover an unbinned example and at the end of the page we include modifications for binned data. This tutorial assumes that you have gone through the standard [likelihood analysis](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_tutorial.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the Data\n",
    "\n",
    "For this thread the original data were extracted from the [LAT data server](http://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi) with the following selections (these selections are similar to those in the paper):\n",
    "\n",
    "```\n",
    "Search Center (RA,Dec) = (251.2054,57.5808)\n",
    "Radius = 30 degrees\n",
    "Start Time (MET) = 322963202 seconds (2011-03-28T00:00:00)\n",
    "Stop Time (MET) = 323568002 seconds (2011-04-04T00:00:00)\n",
    "Minimum Energy = 100 MeV\n",
    "Maximum Energy = 300000 MeV\n",
    "```\n",
    "\n",
    "You will need the following files:\n",
    "```\n",
    "L181102105258F4F0ED2772_PH00.fits\n",
    "L181102105258F4F0ED2772_SC00.fits\n",
    "gll_iem_v07.fits\n",
    "iso_P8R3_SOURCE_V2.txt\n",
    "```\n",
    "\n",
    "Run the code cell below to download them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L181102105258F4F0ED2772_PH00.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/pyLikelihood/L181102105258F4F0ED2772_SC00.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/gll_iem_v07.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/iso_P8R3_SOURCE_V2.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir data\n",
    "!mv *.fits *.txt ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perform Event Selections\n",
    "\n",
    "You could follow the unbinned likelihood tutorial to perform your event selections using **gtlike*, **gtmktime**, etc. directly from the command line, and then use pylikelihood later.\n",
    "\n",
    "But we're going to go ahead and use python. The `gt_apps` module provides methods to call these tools from within python. This'll get us used to using python.\n",
    "\n",
    "So, let's jump into python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gt_apps as my_apps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first run **gtselect** (`filter` in python):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.filter['evclass'] = 128\n",
    "my_apps.filter['evtype'] = 3\n",
    "my_apps.filter['ra'] = 251.2054\n",
    "my_apps.filter['dec'] = 57.5808\n",
    "my_apps.filter['rad'] = 10\n",
    "my_apps.filter['emin'] = 100\n",
    "my_apps.filter['emax'] = 300000\n",
    "my_apps.filter['zmax'] = 90\n",
    "my_apps.filter['tmin'] = 322963202\n",
    "my_apps.filter['tmax'] = 323568002\n",
    "my_apps.filter['infile'] = './data/L181102105258F4F0ED2772_PH00.fits'\n",
    "my_apps.filter['outfile'] = './data/SwiftJ1644_filtered.fits'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.filter.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to find the GTIs. This is accessed within python via the `maketime` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.maketime['scfile'] = './data/L181102105258F4F0ED2772_SC00.fits'\n",
    "my_apps.maketime['filter'] = '(DATA_QUAL>0)&&(LAT_CONFIG==1)'\n",
    "my_apps.maketime['roicut'] = 'no'\n",
    "my_apps.maketime['evfile'] = './data/SwiftJ1644_filtered.fits'\n",
    "my_apps.maketime['outfile'] = './data/SwiftJ1644_filtered_gti.fits'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.maketime.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Livetime Cube and Exposure Map\n",
    "\n",
    "Let's compute the livetime cube and exposure map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Livetime Cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.expCube['evfile'] = './data/SwiftJ1644_filtered_gti.fits'\n",
    "my_apps.expCube['scfile'] = './data/L181102105258F4F0ED2772_SC00.fits'\n",
    "my_apps.expCube['outfile'] = './data/SwiftJ1644_ltCube.fits'\n",
    "my_apps.expCube['zmax'] = 90\n",
    "my_apps.expCube['dcostheta'] = 0.025\n",
    "my_apps.expCube['binsz'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.expCube.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exposure Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from GtApp import GtApp\n",
    "\n",
    "expCubeSun = GtApp('gtltcubesun','Likelihood')\n",
    "expCubeSun.command()\n",
    "my_apps.expMap['evfile'] = './data/SwiftJ1644_filtered_gti.fits'\n",
    "my_apps.expMap['scfile'] ='./data/L181102105258F4F0ED2772_SC00.fits'\n",
    "my_apps.expMap['expcube'] ='./data/SwiftJ1644_ltCube.fits'\n",
    "my_apps.expMap['outfile'] ='./data/SwiftJ1644_expMap.fits'\n",
    "my_apps.expMap['irfs'] = 'CALDB'\n",
    "my_apps.expMap['srcrad'] = 20\n",
    "my_apps.expMap['nlong'] = 120\n",
    "my_apps.expMap['nlat'] = 120\n",
    "my_apps.expMap['nenergies'] = 37"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.expMap.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate XML Model File\n",
    "\n",
    "We need to create an XML file with all of the sources of interest within the Region of Interest (ROI) of SwiftJ1644 so we can correctly model the background.\n",
    "\n",
    "We'll use the user contributed tool `make4FGLxml.py` to create a model file based on the LAT 8-year LAT catalog. You'll need to download the FITS version of this file at http://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/ and get the `make4FGLxml.py` tool from the user-contributed software page and put them both in your working directory.\n",
    "\n",
    "Also make sure you have the most recent galactic diffuse and isotropic model files, which can be found at http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.\n",
    "\n",
    "The catalog and background models are also packaged with your installation of the ScienceTools, which can be found at: `$FERMI_DIR//refdata/fermi/galdiffuse/`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/make4FGLxml.py\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/gll_psc_v19.fit\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/gll_iem_v07.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv *.fit* ./data/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have all of the files we need, we can generate your model file in python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from make4FGLxml import *\n",
    "\n",
    "mymodel = srcList('./data/gll_psc_v19.fit','./data/SwiftJ1644_filtered_gti.fits','./data/SwiftJ1644_model.xml')\n",
    "mymodel.makeModel('./data/gll_iem_v07.fits', 'gll_iem_v07', './data/iso_P8R3_SOURCE_V2.txt', 'iso_P8R3_SOURCE_V2_v1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute the diffuse source responses\n",
    "\n",
    "The [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt) tool will add one column to the event data file for each diffuse source.\n",
    "\n",
    "The diffuse response depends on the instrument response function (IRF), which must be in agreement with the selection of events, i.e. the event class and event type we are using in our analysis.\n",
    "\n",
    "Since we are using SOURCE class, `CALDB` should use the `P8R2_SOURCE_V6` IRF for this tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gt_apps as my_apps\n",
    "\n",
    "my_apps.diffResps['evfile'] = './data/SwiftJ1644_filtered_gti.fits'\n",
    "my_apps.diffResps['scfile'] = './data/L181102105258F4F0ED2772_SC00.fits'\n",
    "my_apps.diffResps['srcmdl'] = './data/SwiftJ1644_model.xml'\n",
    "my_apps.diffResps['irfs'] = 'CALDB'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_apps.diffResps.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the Likelihood Analysis\n",
    "\n",
    "First, import the pyLikelihood module and then the UnbinnedAnalysis functions. For more details on the pyLikelihood module, check out the [pyLikelihood Usage Notes](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_usage_notes.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyLikelihood\n",
    "from UnbinnedAnalysis import *\n",
    "\n",
    "obs = UnbinnedObs('./data/SwiftJ1644_filtered_gti.fits','./data/L181102105258F4F0ED2772_SC00.fits',expMap='./data/SwiftJ1644_expMap.fits',expCube='./data/SwiftJ1644_ltCube.fits',irfs='CALDB')\n",
    "like = UnbinnedAnalysis(obs,'./data/SwiftJ1644_model.xml',optimizer='Minuit')\n",
    "like.tol = 0.1\n",
    "likeobj = pyLike.Minuit(like.logLike)\n",
    "like.fit(verbosity=0,covar=True,optObject=likeobj)\n",
    "likeobj.getQuality()\n",
    "like.logLike.writeXml('./data/Swift_fit1.xml')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This output corresponds to the `MINUIT` fit quality. A \"good\" fit corresponds to a value of `fit quality = 3`; if you get a lower value it is likely that there is a problem with the error matrix. Now we try with NewMinuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "like2 = UnbinnedAnalysis(obs,'./data/Swift_fit1.xml',optimizer='NewMinuit')\n",
    "like2.tol = 0.0001\n",
    "like2obj = pyLike.NewMinuit(like2.logLike)\n",
    "like2.fit(verbosity=0,covar=True,optObject=like2obj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(like2obj.getRetCode())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you get anything other than `0` here, then NewMinuit didn't converge.\n",
    "\n",
    "We can start by deleting sources with low or negative TS, which tend to hinder convergence. First, we delete sources with TS levels below 10 and run the fit again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sourceDetails = {}\n",
    "\n",
    "for source in like2.sourceNames():\n",
    "    sourceDetails[source] = like.Ts(source)\n",
    "\n",
    "for source,TS in sourceDetails.iteritems():\n",
    "    print source, TS\n",
    "    if (TS < 10):\n",
    "        print(\"Deleting...\")\n",
    "        like2.deleteSource(source)\n",
    "\n",
    "like2.fit(verbosity=0,covar=True,optObject=like2obj)\n",
    "\n",
    "print(like2obj.getRetCode())\n",
    "\n",
    "like2.logLike.writeXml('./data/Swift_ts10.xml')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have achieved convergence, we need to add the SwiftJ1644 source to the top of `Swift_ts10.xml` model file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```xml\n",
    "<source name=\"SwiftJ1644\" type=\"PointSource\">\n",
    " <spectrum type=\"PowerLaw2\">\n",
    "  <parameter free=\"true\" max=\"10000.0\" min=\"0.0001\" name=\"Integral\" scale=\"1e-07\" value=\"1.0\"/>\n",
    "  <parameter free=\"true\" max=\"5.0\" min=\"0.0\" name=\"Index\" scale=\"-1.0\" value=\"2.0\"/>\n",
    "  <parameter free=\"false\" max=\"500000.0\" min=\"20.0\" name=\"LowerLimit\" scale=\"1.0\" value=\"100.0\"/>\n",
    "  <parameter free=\"false\" max=\"500000.0\" min=\"20.0\" name=\"UpperLimit\" scale=\"1.0\" value=\"300000.0\"/>\n",
    " </spectrum>\n",
    " <spatialModel type=\"SkyDirFunction\">\n",
    "  <parameter free=\"false\" max=\"360.0\" min=\"-360.0\" name=\"RA\" scale=\"1.0\" value=\"251.2054\"/>\n",
    "  <parameter free=\"false\" max=\"90.0\" min=\"-90.0\" name=\"DEC\" scale=\"1.0\" value=\"57.5808\"/>\n",
    " </spatialModel>\n",
    "</source>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the Swift source in the XML file, we can now calculate the upper limit (the paper used an upper energy limit of 10GeV so that's what we are using here)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "help(UnbinnedAnalysis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "like3 = UnbinnedAnalysis(obs,'./data/Swift_ts10TEST.xml',optimizer='Minuit')\n",
    "\n",
    "from UpperLimits import UpperLimits\n",
    "\n",
    "ul = UpperLimits(like3)\n",
    "ul['SwiftJ1644'].compute(emin=100,emax=10000)\n",
    "print ul['SwiftJ1644'].results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this is in ph/cm^2/s and not ergs/cm^2/s."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Binned Data Upper Limits\n",
    "\n",
    "In the case of binned data, one needs to follow the standard [likelihood analysis](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html) or the [python version](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/summed_tutorial.html) in order to generate livetime cubes, counts cubes, source maps and exposure maps.\n",
    "\n",
    "The upper limits steps are nearly identical to the previous section, but one needs to import BinnedAnalysis and use BinnedObs with the proper file format instead:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "import pyLikelihood\n",
    "from BinnedAnalysis import *\n",
    "\n",
    "obs = BinnedObs(srcMaps='file_name',binnedExpMap='file_name',\n",
    "expCube ='file_name',irfs='CALDB')\n",
    "like = BinnedAnalysis(obs,'XML_file_name',optimizer='NewMinuit')\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
